%% Maximize the Likelihood function


clear
delete('Ryan');
diary('Ryan');

cd ../china/data
load hs10_indlist.raw;
indN=length(hs10_indlist);

for w=1:50

folder=sprintf('../china/estimation/hs10/regular/raw_%d',hs10_indlist(w,1))


Name='Ryan is Great';

tic
addpath('/apps/tomlab');
startup

cd (folder)


load data2_q_regular.mat


nExp = X;
nImp = M;
S=N*nExp*nExp;
delta=0.975;

failnum=0;


% tom objects include beta_p, beta_x and EV
beta_p=tom('beta_p',1,1);
beta_x=tom('beta_x',1,1);
beta_c=tom('beta_c',1,1);
xi_tilde=tom('xi_tilde',1,1);
EV=tom('EV',S,1);

U=beta_p*ep + beta_x*switched + beta_c*city_switched + xi_tilde*lambda;

% EV is a tomSym object, so cannot write:
% for i=1:S CSVF(i,1) = U(i,1) + delta*EV(i,1)

CSVF= U + delta*EV;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P



 P = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;



% New Constraint Definition:

W=tomArray(CSVF,[N*X X]);
W=repmat(W,N,1);
W=W(:);
W=repmat(W,X,1);
W=tomArray(W,[N X N*X*X]);
Y=  log( sum( exp( W),2) );
TransProbs=reshape(TransProbs,N,N*X*X);
Z=sum( Y.*TransProbs,1);
Z=Z(:);




objective = sum ( log ( P ( impchoicestate) ) );

constraints = {EV==Z};


EV_guess=ones(S,1);

beta_x_guess=-1;
beta_c_guess=-1;
guess = struct('beta_p',0,'beta_x',beta_x_guess,'beta_c',beta_c_guess,'xi_tilde',0,'EV',EV_guess);
options = struct;
options.name = Name;
options.prilev=3;
options.use_H=false;
options.use_d2c=false;

Prob = sym2prob('minlp',-objective,constraints,guess,options);

%Prob.DUNDEE.optPar(20) = 1;
%Prob.optParam.IterPrint = 1;
Prob.KNITRO.options.Feastol_Abs=0.0001;
Prob.KNITRO.options.Opttol_Abs=0.0001;
Prob.PriLevOpt=3;

Result = tomRun('knitro',Prob,2);
toc

b=Result.x_k;
beta_c_sol=b(S+1);
beta_p_sol=b(S+2);
beta_x_sol=b(S+3);
xi_tilde_sol=b(S+4);
EV_sol=b(1:S);

[beta_p_sol beta_x_sol beta_c_sol xi_tilde_sol]

save ('beta_regular', 'beta_p_sol','beta_x_sol','beta_c_sol','xi_tilde_sol'); 
save ('EV','EV_sol');

w

cd ../china/estimation/hs10

end

diary off